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We present a detailed numerical study of ground state and finite temperature spin and charge 
dynamics of the two-dimensional Kondo lattice model with hopping t and exchange J. Our nu- 
merical results stem from auxiliary field quantum Monte Carlo simulations formulated in such a 
way that the sign problem is absent at half-band filling thus allowing us to reach lattice sizes up 
to 12 x 12. At T = and antiferromagnetic couplings, J > 0, the competition between the RKKY 
interaction and Kondo effect triggers a quantum phase transition between antiferromagnetically or- 
dered and magnetically disordered insulators: J c /t = 1.45 ± 0.05. At J < the system remains 
an antiferromagnetically ordered insulator and irrespective of the sign of J, the quasiparticle gap 
scales as |J|. The dynamical spin structure factor, S(q,uj), evolves smoothly from its strong cou- 
pling form with spin gap at q = (tv, it) to a spin wave form. For J > 0, the single particle spectral 
function, A(k,co), shows a dispersion relation following that of hybridized bands as obtained in the 
non-interacting periodic Anderson model. In the ordered phase this feature is supplemented by 
shadows thus allowing an interpretation in terms of coexistence of Kondo screening and magnetic 
ordering. In contrast, at J < the single particle dispersion relation follows that of non-interacting 
electrons in a staggered external magnetic field. At finite temperatures spin, Ts, and charge, Tc, 
scales are defined by locating the maximum in the charge and spin uniform susceptibilities. For 
weak to intermediate couplings, Ts marks the onset of antiferromagnetic fluctuations - as observed 
by a growth of the staggered spin susceptibility- and follows a J 2 law. At strong couplings Ts scales 
as J. On the other hand Tc scales as J both in the weak and strong coupling regime. At and 
slightly below Tc we observe i) the onset of screening of the magnetic impurities, ii) a rise in the 
resistivity as a function of decreasing temperature, iii) a dip in the integrated density of states at 
the Fermi energy and finally iv) the occurrence of hybridized bands in A(k,u>). It is shown that in 
the weak coupling limit, the charge gap of order J is formed only at Ts and is hence of magnetic 
origin. The specific heat shows a two peak structure. The low temperature peak follows Ts and is 
hence of magnetic origin. Our results are compared to various mean-field theories. 
PACS numbers: 71.27.+a, 71.10.-w, 71.10.Fd 



I. INTRODUCTION 

The Kondo lattice model (KLM) as well as the periodic Anderson model (PAM) are the prototype Hamiltonians 
to describe heavy fermion materials jjj and Kondo insulators ||. The physics under investigation is that of a lattice 
of magnetic impurities embedded in a metallic host. The symmetric PAM reads: 

h pam = y: eCk)ct a%rT -vj: (c.j-. + 4 Cla ) + Uf x: <v T - 1/2) (4 _ i /2 ) . w 

k.tj i,cr i 

The unit cell, denoted by i, contains an extended and a localized orbital. The fcrmionic operators ct (ft ) create 

electrons on extended (localized) orbitals with wave vector k and z— component of spin a. The overlap between 
extended orbitals generates a conduction band with dispersion relation e(k). There is a hybridization matrix element, 
V , between both orbitals in the unit-cell and the Coulomb repulsion- modeled by a Hubbard U /- is taken into account 
on the localized orbitals. In the limit of large Uf, charge fluctuations on the localized orbitals are suppressed and the 
PAM maps onto the KLM |§: 

H klm = E e(k)cl a%a + JJ2SS-S(. (2) 

k.o i 

Here SS = h ; d ct, ,/C? where a are the Pauli s — 1/2 matrices. A similar definition holds for SI. A magnetic 
energy scale J = 8V 2 /U emerges and there is a constraint of one electron per localized orbital. Although this constraint 
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forbids charge fluctuations on the localized orbitals, those fluctuations are implicitly taken into account leading to 
the above form and sign of the exchange interaction. On the other hand, when charge fluctuations on the localized 
orbitals are absent, the exchange interaction follows from Hund's rule and is ferromagnetic. The ferromagnetic KLM 
has attracted much attention in conjunction with manganites [|J. In this article we will consider both ferromagnetic 
and antiferromagnetic exchange interactions with emphasis on the antiferromagnetic case. 

The physics of the single impurity Anderson and Kondo models at J/t > is well understood || . In the temperature 
range J < T < U charge is localized on the /-orbital, but the spin degrees of freedom are essentially free thus leading 
to a Curie- Weiss law for the impurity spin susceptibility. Below the Kondo temperature Tk oc e /e~ ly IJN ( £ f"> the 
impurity spin is screened by the conduction electrons. Here, Sf is the Fermi energy and N(ef) the density of states 
taken at the Fermi energy. The transition from high to low temperatures is non-perturbative and corresponds to the 
Kondo problem with the known resistivity minimum |jj and orthogonality catastrophe . At low temperatures Tr- 
is the only energy scale in the problem. 

A lattice of magnetic impurities introduces new energy scales. In the spin sector, the Ruderman-Kittel-Kasuya- 
Yosida (RKKY) interaction || couples impurity spins via polarization of the conduction electrons. This interaction 
takes the form of a Heisenberg model with exchange J e ff(q) oc — J 2 Kex(q, u> = 0) where x(?i a -') corresponds to the 
spin susceptibility of the conduction electrons. Since this interaction favors magnetic ordering, it freezes the impurity 
spins and hence competes with the Kondo effect. By comparing energy scales one expects the RKKY interaction 
(Kondo effect) to dominate at weak (strong) couplings. As suggested by Doniach |J, this leads to a quantum phase 
transition between ordered and disordered magnetic phases. 

As a function of dimension, contrasting results are obtained for the PAM and KLM. We first consider the limit 
of large dimensions [ |To| , |TT[ ] and the Gutzwiller approximation JT^ ] . The Gutzwiller approximation leads to an non- 
interacting PAM with renormalized hybridization V. At half-filling an insulating state is obtained, with quasiparticle 
gap ~ e - 1 / 2 J N ( e f) i n the large Uf limit. Both the Gutzwiller and dynamical mean-field approaches yield charge and 
spin gaps equal to each other. As a function of temperature, optical and quasiparticle gaps start appearing at an 
energy scale ~ e - 1 / 2JAr ( £ /) [fij|. I n the doped phase, the Luttinger volume includes the /-electrons, and due to the 
renormalization of the hybridization, the effective mass of charge carriers is enhanced. The above quoted results stem 
from calculations for the PAM. However, similar results are obtained in the framework of the KLM at J/t << I in 
the limit of large dimensions Jff| . The above approximations predict an instability to magnetic ordering in the large 
Uf or small J limit. The occurrence of this instability has been observed in the framework of quantum Monte Carlo 
(QMC) simulations of the PAM in two dimensions 

HQ. In the one-dimensional good understanding of the 

phase diagram of the KLM as a function of electronic density and coupling has been achieved Jll],[l6| . In particular 
at half-filling, a spin liquid phase is obtained irrespective of the value of J/t. In the weak coupling limit the spin gap 
follows a Kondo form, whereas the charge gap tracks J. 

In this article we present a detailed numerical study of ground-state and finite-temperature properties of the half- 
filled KLM in intermediate dimensions, d = 2. Our T = simulations arc aimed at understanding the competition 
and interplay of the Kondo effect and RKKY interaction. Our finite temperature simulations provide insight into the 
temperature evolution of spin and charge degrees of freedom. 

Our main results and structure of the article is as follows. Details of the numerical technique are presented in the 
next section. We use a path integral auxiliary field quantum Monte-Carlo (QMC) method Our approach is based 
on a simple technical innovation which allows to avoid the sign problem at least at half-band filling where the model 
is particle-hole symmetric. Both finite and zero temperature versions of the algorithm are presented. In both cases 
imaginary time displaced correlations functions can be computed. The continuation to real time is then carried out 
via the use of the maximum entropy (ME) method [fill. We note that the algorithms may be applied irrespective of 
the sign of J. 

In section III ground state equal time and dynamical properties of the ferromagnetic and antiferromagnetic KLM 
are presented. Our main results include the following, i) In the spin sector, a quantum phase transition between 
antiferromagnetically ordered and disordered ground states occurs at J/t = 1.45±0.05. The dynamical spin structure 
factor is analyzed across the transition. As a function of decreasing values of J/t, the spin gap at the antiferromagnetic 
wave vector closes and the magnon spectrum evolves towards a spin-wave form. This spin wave form persists for 
ferromagnetic couplings since in the limit J/t — > oo the model maps onto the s = 1 antiferromagnetic Heisenberg 
model. Our results at J/t > are compared to a bond-operator mean field theory of the Kondo necklace model, ii) 
In the charge sector, the system remains an insulator. To a first approximation, the quasiparticle gap tracks J both 
in the antiferromagnetic and ferromagnetic KLM. For all values of J/t > the single-particle spectral function shows 
a feature whose dispersion relation follows the one obtained in the non-interacting PAM. In a mean-field approach, 
this feature results solely from Kondo screening of the magnetic impurities. In the magnetically ordered phase, this 
feature is supplemented by shadow bands. Thus and as confirmed by a mean-field approach, the spectral function in 
the ordered phase may only be understood in terms of the coexistence of Kondo screening and the RKKY interaction. 
On the other hand, at J/t < where Kondo screening is absent the single particle dispersion relation follows that of 
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free electrons in a external staggered magnetic field. 

Section |[y]is devoted to finite temperature properties of the KLM. We define charge, Tc, as well as spin, Tg, scales 
from the location of the maximum in the charge and spin susceptibilities. In the weak and strong coupling limit, the 
charge scale tracks J. On the other hand the spin scale - as expected form the energy scale associated to the RKKY 
interaction - follows a J 2 law up to intermediate couplings. At strong couplings Tg cx J. Since Tc corresponds to 
the energy scale at which a minimum in the resistivity is observed, we conclude that it describes the energy scale at 
which scattering is enhanced due to the screening of the impurity spins. Furthermore a reduction of the integrated 
density of states at the Fermi level is observed at Tc- The spin scale up to intermediate couplings (i.e. J/t < W where 
W corresponds to the bandwidth) marks the onset of short-range antiferromagnetic correlations. This is confirmed 
by the calculation of the staggered spin susceptibility which shows a strong increase at Tg. In the weak coupling 
limit, it is shown that the quasiparticle gap of magnitude oc J is formed only at the magnetic energy scale Tg and 
is thus of magnetic origin. In the temperature range Tg < T < Tc hybridized band are seen in the single particle 
spectral function with quasiparticle gap lying beyond our resolution. Finally, the specific heat is computed and shows 
a two-peak structure. The low energy peak tracks the spin scale and is hence of magnetic origin. 

In the last section, we discuss our results as well as links with experiments. 



II. AUXILIARY FIELD QUANTUM MONTE-CARLO ALGORITHM FOR THE KONDO LATTICE 

MODEL 

Auxiliary held QMC simulations of the KLM as well as the two-impurity Kondo model have already been carried 
out by Fye and Scalapino as well as by Hirsch and Fye However, their formulation leads to a sign problem 

even in the half-filled case where the model is invariant under a particle-hole transformation. In this section we present 
an alternative formulation of the problem which is free of the sign problem in the particle-hole symmetric case. In 
order to achieve our goal, we take a detour and consider the Hamiltonian: 



-l 2 

< J i.cr i.(7 1,0 i.t 



(3) 



As we will see below, at vanishing chemical potential this Hamiltonian has all the properties required to formulate a 
sign-free auxiliary held QMC algorithm. Here, we are interested in ground-state properties of H which we obtain by 
filtering out the ground state \Hfo) by propagating a trial wave function \^t) along the imaginary time axis: 

(tt |O|tt ) (y T \e-° H Oe- eH \y T ) 
(*o|*o> (* T |e- 20ff |* T ) 1 ' 

The above equation is valid provided that (^t\^o) 7^ and O denotes an arbitrary observable. 
To see how H relates to Hklm we compute the square in Eq. (^) to obtain: 

H = Ye(k)c\ c r +jYsS-Sl-^-Y(ci cl h n h n + H.C.) + ~ Y ( r£nl - n% - n{ ). (5) 

^-^/ y ' k,<7 k,rr Z-^i i i 4 Z-^i \ i,a- i,—tr Jt ,-<7 J h<7 J 4 \ 1 1 1 i J w 

As apparent, there are only pair-hopping processes between the /- and c-sites. Thus the total number of doubly 
occupied and empty /-sites is a conserved quantity: 

[ ff > £(* - n { T )( 1 - ™L) + n ?, T n L] = °- (6) 

i 

If we denote by Q n the projection onto the Hilbert space with YllO- ~ n ^ t )(l — n ^ ) + nl nl = n then: 

JN 

HQ = Hklm + — (7) 

since in the Qo subspace the /-sites are singly occupied and hence the pair-hopping term vanishes. Thus, it suffices 
to choose 

Qo|*t) = |*t) (8) 
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to ensure that 

(y T \e~ eH Oe~ eH \->Z T ) _ (q T \e- eHK ^'Oe- eH »-^<\y T ) 



(9) 



It is interesting to note that there is an alternative route to obtain the KLM. Instead of projecting onto the Q 
Hilbcrt space, we can project onto the Qn Hubert space by suitably choosing the trial wave function. 

HQ N = £ e{k)c\ A - £ £ (iJ,-A-«fl° + H ' c -) +lE ( n H ~ n l~ n ■ ( 10 ) 

k.a i.<7 i 

Since in the Qn subspace the /-sites are doubly occupied or empty, the exchange term SS ■ Si vanishes. To see the 
relation with the KLM, we define the spin-1/2 operators: 

which operate on the states: | f = — (— l) lx+Iy /J t /J |0) and | t — |0) as well as the fermion operators: 



= d., 4 = (-l)*-+*»^; r (12) 



HQn = ^4,^ + f E Or ^ + W) + ^ J + ^ d3) 



With those definitions, 
which is nothing but the KLM. 

A. The basic formalism 

Having shown the relationship between H and Hklm we now discuss some technical aspects of the QMC evaluation 
of {^ T \e~ @H Oe~ @H \^ T )/{^T\e~ 2@H \^T)- With the use of the Trotter formula we obtain: 

M 

(* T |e- 26H |* T > = (* r | Yl e ArHt e- ATHj \y T ) + 0(At 2 ) (14) 

T=l 

Here H t — — o- * C J a + H.c, Hj = — j^2^S^ ■ Si, and MAt = 20. Strictly speaking, the systematic error 

produced by the above Trotter decomposition should be of order At. However, if the trial wave function as well as 
H t and Hj are simultaneously real representable, it can be shown that the prefactor of the linear At error vanishes 

Since we will ultimately want to integrate out the fcrmionic degrees of freedom, we carry out a Hubbard- 
Stratonovitch (HS) decomposition of the perfect square term p3[ : 

e -ATH J= -Q e A ^/ 4 (E,i.4.+H.c.) 2 

i 

= n( E l(l)e VX7JrMl) ^ c ^ +H - c -+0(At 4 )] , 

r y=±i,±2 / 

where the fields rj and 7 take the values: 

7(±1) = 1 + V6/3, 7 (±2) = 1 - n/6/3 



(15) 



T)(±l) =±4/2(3-^6), r?(±2) = ±W2 (3 + V6 
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As indicated, this transformation is approximate and produces on each time slice a systematic error proportional to 
At 4 . This amounts to a net systematic error of order MAr 4 ~ 28Ar 3 which for constant values of the projection 
parameter is an order smaller that the error produced by the Trotter decomposition. 

The trial wave function is required to be a Slater determinant factorizable in the spin indices: 

|* r ) = |4) ® |**) with = (E<* P «,J 10). (16) 

y=l \ x / 

Here, we have introduced the notation x = (i.n) where i denotes the unit cell and n the orbital (i.e. at- , = cl 

and at- , = ft ). It is convenient to generate I^S,) from a single particle Hamiltonian HZ = at (/i?) a,, 

which has the trial wave function as non-degenerate ground state. To obtain a trial wave function which satisfies the 
requirements QoI^t) = \^t) we are forced to choose Ho of the form: 

H »= E M,. c - + H - c -) + ? e ^' ? (4t^t ~ ( 17 ) 

which generates a Neel state (Q = (7r, 7r)) on the localized orbitals. To obtain a non-degenerate ground state, we 
impose the dimerization 

with 5 « t. 

We are now in a position to integrate out the fermionic degrees of freedom to obtain: 

(* T |e- 2e "|*r> = £ (l[ 7 (l lT )) IJdet (i** e -±rfJ ( r )pa \ (ig) 

« \i,T / o- V r=l / 

where the matrices T and J(t) are defined via: 

H t = E £ (^4, CT c fe^ = E 4, A»ay,» 

x,y,a 

The HS field / has acquired a space, i, and time, r, index. 

The basic ingredients to compute observables are equal-time Green functions. They are given by: 

M/2 M/2+1 

£/> = J] e~ ATi e 3 ^P° U< = \{ e- AT *e J V, and 

t=1 r=M 



(n?, T 7fe T ))n CT det (t/<EP 



Pr(Z) = — ^-r V j -At - (21) 

£{/} (n?, T 7(^))n CT det ([/<,{/>) 

Since, for a given set of HS fields, we are solving a free electron problem interacting with an external field a Wick 
theorem applies. Hence from the knowledge of the the single particle Green function at fixed HS configuration we 
may evaluate all observables. Imgaginary time displaced correlation functions may equally be calculated E4J23]. 



5 



We are left with the summation over the HS fields which we will carry out with Monte-Carlo methods. In order to 
do so without further complication, we have to be able to interpret Pr(Z) as a probability distribution. This is possible 
only provided that Pr(7) > for all HS configurations. In the particle-hole symmetric case the above statement is 
valid. Starting from the identity: 

Tr (e-^o H M =1 e- A - H ie H l^) 
det (UftUfA = lim — ^ ] —— . >- (22) 

we can carry out a particle-hole transformation: 

4, T "> and -> (23) 

Here, H% = J2x y a x,a- r ^x,ya.y,a and Hj{t) = y a x,<jj{ T )x,y a y-iy Since Eq. ( p3|) corresponds to a canonical transfor- 
mation, the trace remains invariant and Hq, H\ as well as Hj(t) map onto H^, and Hj(t) respectively. Thus we 
have shown that: det 

( u u u u) = det ( u u u u) which leads to Pr ^) - for a11 values 01 the HS nelds - Awa y from 

half- filling (which would correspond to adding a chemical potential term in Hq), particle hole-symmetry is broken and 
Pr(^) may become negative. This leads to the well known sign-problem R 

For the Monte-Carlo sampling of the probability distribution Pr(l), we adopt a sequential single spin-flip algorithm. 
The details of the upgrading procedure as well as of the numerical stabilization of the code are similar to those used 
for auxiliary field QMC simulations of the Hubbard model [Eq] . 



B. Optimizing the algorithm 

The above straightforward approach for the QMC simulation of H turns out to be numerically inefficient. The 
major reason for this stems from the choice of the trial wave function. The coupled constraints i) Qo\^t) = \^t) and 
ii) I^&t) is a Slater determinant factorizable in the spin indices make it impossible to choose a spin-singlet trial wave 
function (the trial wave function generated by the single particle Hamiltonian Hq of Eq. ( |l7|) orders the /—electrons 
in a Neel states which is not a spin singlet). Since we know that the ground state of the KLM on a finite-size system is 
a spin singlet (|^J2^] , we have to filter out all the spin excited states from the trial wave function to obtain the ground 
state. This is certainly not a problem when we are investigating the physics of a problem with a large spin-gap as is 
the case in the limit J/t >> 1. However, in the limit of small J/t long-range magnetic order is present and hence one 
expects finite-size spin-gap to scale as v s /L where v s is the spin velocity and L the linear size of the system. In this 
case, starting with a spin-singlet trial wave function is important to obtain reliable convergence p4j| . 

In order to circumvent the above problem, we relax the constraint Qo\^t) — \^t) and add a Hubbard term for the 
/-sites to the Hamiltonian. 



fc,cr i 



yd h + ft c~ 



tf/EWt-VaXn^-l^). ( 24 ) 



This Hamiltonian is again block diagonal in the Q n subspaces. During the imaginary time propagation, the components 
Qnl^r) of the trial wave function will be suppressed by a factor e~ eUfn ^ 2 in comparison to the component QoI^t). 
The usual procedure to incorporate the Hubbard term in the QMC simulation relies on Hirsch's HS transformation 




H-^EKi-dK-^JI (25) 

E ex p 5 E s ?( n ?,T _n a) • 

S1,...,SJV=±1 \ ? / 



1 It is clear that by choosing Hq = Hq thus leading to = P^ would produce positive values of Pr(£) for all HS configurations 
and irrespective of particle-hole symmetry. This stands in analogy to the absence of sign-problem in the attractive Hubbard 
model. However, this choice of the trial wave function is incompatible with the requirement Qo\^t) = |^t). 



G 



where cosh(a) = exp (AtII/2). As apparent from the above equation, for a fixed set of HS fields, s% . 
symmetry is broken. Clearly SU(2) spin symmetry is restored after summation over the HS fields 
Alternatively, one may consider [^9| 



s N , SE/(2)-spin 



exp -AtU^2 



t 



(26) 



C 



E 



exp za 



,SN=±1 



»,1 



where cos(a) = exp(-Art//2) and C = exp(ArUN/4) /2 N . With this choice of the HS transformation 5/7(2) spin 
invariance is retained for any given HS configuration. Even taking into account the overhead of working with complex 
numbers, one of the authors has argued |]30| that this choice of HS transformation produces a more efficient code. 

Having relaxed the condition Qo\^t) = \^t) we are now free to choose a spin singlet trial wave function which we 
generate from: 



F = ]T e (fc)4 



J 

4 ^ — ' * 

i,<j 



V(ct h + f-l c- ) 



(27) 



which is nothing but the non-interacting PAM with hybridization V = J/4. The ground state at half-filling is clearly 
a spin singlet. With this choice of the trial wave function, and the Hubbard-Stratonovitch transformation of Eq. ( |2^ ) 

the particle-hole transformation of Eq. (23) maps det ff/jr 
occurs at half-filling. 



on det (u^i s U^i s ^j- Hence, no sign problem 



L = A, J/t = 1.6 



A: S xy (Q), 5 2 |* r ) ^0 
.: S(Q), 5 2 |* r ) = 
O: S Z (Q), 5 2 |* r ) ^0 



Q I I * <S I I I I I I 

5 10 15 20 25 30 

2Qt 

FIG. 1. Spin-spin correlations as a function of the projection parameter O. Here, S(Q) = %{S*(Q) ■ S^(~Q)), 
Sl(Q) = 4{5|(Q) • S£(-Q)), and S f xy (Q) = 2 ((S[(Q) ■ S[{-Q)) + 0f(Q) ■ S^(-Q))). The trial wave function with 

S^I^t) 7^ ( S 2 |*t) = 0) corresponds to the ground state of the Hamiltonian in Eq. (gTj) (Eq. (0)). In the large Q 
limit, the results are independent on the choice of the trial wave function. In particular, starting from a broken symmetry state 
the symmetry is restored at large values of Ot. For this system, the spin gap is given by A sp = 0.169 ± 0.004 J31J. Starting 
with a trial wave function with S 2 |^1/t) 7^ 0, convergence to the ground state follows approximatively the form: a + be~ A ' p20 . 
The solid lines correspond to a least square fit to this form. 

Fig. [l] demonstrates the importance of using a spin singlet trial wave function. Starting from a Neel order for the 
f-electrons, convergence to the ground state follows approximatively e _As *> 2e> where A sp corresponds to the spin-gap. 
When the spin gap is small, convergence is poor and the remedy is to consider a spin singlet trial wave function. 

Having optimized the trial wave function we now consider convergence as a function of J/t. As apparent from Fig. |2| 
for small values of J/t increasingly large projection parameters are required to obtain convergence. The origin of this 
behavior may be traced back to the energy scale of the RKKY interaction which follows a J 2 law. At J/t = 0.4, 
2<dt ~ 40 is enough to obtain convergence whereas at J/t = 0.2, a value of 2<dt ~ 170 is required. 
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S ff (Q = (7t,tt))/L 2 , (n) = 2, L = 6 




O: J/t 


= 0.2 


A: J/t 


= 0.4 


v- J/t 


= 0.8 


*:J/t 


= 1.2 


.: J/t 


= 1.4 


0.08 


0.09 



0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 

1/(26*) 

FIG. 2. Spin structure factor at Q = (7T, ?r) for the /-electrons {S ff (Q)) at various values of J/t and as a function of the 
projection parameter Qt. Here we consider a spin singlet trial wave function. 

The systematic error produced by the Trotter decomposition scales as (At) 2 . This behavior is shown in Fig. (||). 
All our calculation were carried out at values of At small enough so as to neglect this systematic error. 

(n) = 2, L = 4, J/t = 1.6, 6* = 11 



0.22 
0.20 
0.18 - 
0.16 
0.14 
0.12 
0.10 - 



-@ — @- 



O: Sff(Q)/L 2 

Q = (tt.tt) 
A: 5(Q)/L 2 



0.3 

At* 



— A~ 
i 

0.5 



FIG. 3. Systematic error produce by the Trotter decomposition. In our simulations, we have used At — 0.1 and At — 0.2. 
Here, S(Q) corresponds to the spin structure factor of the total spin at Q = (it, 7r). 



C. Ferromagnetic exchange 



Until now, we have implicitly considered an antiferromagnetic exchange, J > 0. It is straightforward to generalize 
the above case to a ferromagnetic one. The only point to take care of is the choice of the trial wave function in order 
to avoid the sign problem. In this case the non-interacting Hamiltonian which generates the trial wave function has 
to be invariant under the particle-hole transformation: 



(-1) 



i x +i. 



z tl and /J - +(-i)*.+H,/ ?iJ 



(28) 



Note that in comparison to Eq. fl23| ) there is an overall sign difference in the particle-hole transformation of the 

f-operators. With this condition one has: det [U^i s U^i = det S U^ so that no sign problem occurs. The 
trial wave function is thus generated from the non-interacting Hamiltonian: 



H = V e(k)cl cr - - V (ct /■» +/] ). 



(29) 
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D. Finite temperature algorithm 



The QMC method presented above may be generalized to finite temperatures to compute expectation values of 
observables in the grand-canonical ensemble: 

, , Tr(e-^O) , 

(o) = -^my ( 3 °) 

Since the step from the T — approach to the finite-T algorithm is similar to the one for the standard Hubbard model, 
we refer the reader to the Ref. p6| ]. We note however, that at finite temperatures, the projection onto the Qq subspace 
may only be achieved via the inclusion of the Hubbard term Uf ^^(nl — \/2){nl — 1/2) in the Hamiltonian. At 

this point, it is very convenient to choose the 5'C/(2)-invariant HS decomposition of Eq. ( |26| ) since one can take the 
limit Uf — ► oo by setting a = tt/2. Hence irrespective of the considered temperature, we are guaranteed to be in the 
correct Hilbert space. 



III. SPIN AND CHARGE DEGREES OF FREEDOM AT T = 

The different phases occurring at half-filling are summarized in Fig. ^. All quantities have been extrapolated to the 
thermodynamic limit M]. We have considered sizes ranging from 4 x 4 to 12 x 12 with periodic boundary conditions. 
The staggered moment: 

m s = lim JUs(Q)-S(-Q)) (31) 

indicates the presence of long-range magnetic order. Here, S(Q) — X)j e where S(j) = S*(j) + S c (j) is the 

total spin, Q = (tt, tt) the antiferromagnetic wave vector and L corresponds to the linear size of the system. This 
quantity is maximal at J/t — —oo and vanishes at J c /t ~ 1.45 thus signaling a phase transition. The onset of a spin 
gap, 

A sp = lim E^{S = 1, N p = 2N) -E^{S = 0, N p = 2N), (32) 

L — >oo 

is observed when magnetic order disappears. Here, Eq (S, N p ) is the ground state energy on a square lattice with 
N = L 2 unit cells, N p electrons and spin S. Finally, the system remains an insulator for all considered coupling 
constants. This is supported by a non-vanishing quasiparticle gap, 

A qp = lim E^{S = 1/2, N p = 2N + 1) -E%(S = 0,N p = 2N). (33) 

L — >oo 

We will first discuss the spin degrees of freedom and then turn our attention to charge degrees of freedom. 
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FIG. 4. (a) Staggered moment, m s , spin gap A sp and quasiparticle gap for the ferromagnetic and antiferromagnetic KLM. All 
quantities have been extrapolated to the thermodynamic limit based on results on lattice sizes up to 12 x 12. The data for J > 
stems from Ref. p3|. The staggered moment corresponds to that of the total spin (see Eq. (J3l|) ) . The solid line corresponds 
to the value of the staggered moment for the s — 1 antiferromagnetic model as obtained in a spin wave approximation M . (b) 
Staggered moment of the /— and c— electrons after extrapolation to the thermodynamic limit. 



A. Spin degrees of freedom 



To investigate the spin degrees of freedom we compute the dynamical spin susceptibility, 

S(q,u>) = nJ2 \{n\S(q)\0)\ 2 t(" ~ ~ Eq)). 



(34) 



where the sum runs over a complete set of eigenstates and |0) corresponds to the ground state. This quantity is 
related to the imaginary time spin-spin correlations which we compute with the QMC method [plj : 



(0|5(<f, r) • S(-q)\0) = - I do;e-™S(g». 

7T 



(35) 



Here, S(q,r) = e rH S(q)e~ rH . We use the Maximum Entropy (ME) method to accomplish the above numerically ill 
defined inverse Laplace transform |[8| . 

In the strong coupling limit J — > oo, the model becomes trivial, since each /—spin captures a conduction electron 
to form a singlet. In this limit, the ground state corresponds to a direct product of singlets on the f-c bonds of a 
unit cell. Starting from this state, one may create a magnon excitation by breaking a singlet to form a triplet. In 
second-order perturbation in i/J, this magnon acquires a dispersion relation given by: 



E sp (q) = J — 



16t 2 
~3J~ 



At 2 



(36) 



where "f{q) — cos(q x ) + cos(q y ) |l5| ]. At Q = (7r,7r), E sp (q) is minimal and is nothing but the spin gap. In Fig. |^a, 
we plot the dynamical spin structure factor for J/t = 2.0. The solid bars in the plot correspond to a fit to the above 
strong coupling functional form: a + bj(q). As apparent, this functional form reproduces well the QMC data. We 
note that this magnon mode lies below the particle hole continuum located at 2A qp (see Fig. ||). 
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FIG. 5. Dynamical spin structure factor at T = for the ferromagnetic and antiferromagnetic KLM. We have normalized 
the peak heights to unity. The numbers on the left hand side of the figures correspond to the normalization factor. The vertical 
bars are fits to the data, see text. 

As we approach the antiferromagnetically ordered phase, one expects that the above magnon mode evolves towards 
a spin- wave form: 



£-® = -fVl- 7(5)74. (37) 

As apparent from Fig. ||d as one approaches J c the spin gap vanishes and the magnon mode softens around q = 0. In 
the antiferromagnetic phase (see Fig. ||c) the data follow well the above spin-wave form. 

In the limit of large ferromagnetic couplings, the model maps onto the S = 1 antiferromagnetic Heisenberg model. 
At J/t = —00, the ground state is macroscopically degenerate since the f-c bonds are effectively decoupled and 
occupied by a triplet with arbitrary z-component of spin. This degeneracy is lifted in second-order perturbation 
theory, yielding aS = l antiferromagnetic Heisenberg model: 

0+2 

He ff = —Y. s rSr (38) 

ill) 

Here, S, = E^iA^' *fo = 4t4< = A (4 T 4 + 4i4)' and ^ = 4A COTres P° nd 
to the s — 1 Pauli spin matrices. The magnetic excitations are clearly spin waves as confirmed by the QMC data 
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of Fig. |5Jd. In the limit of large negative J, the staggered moment should scale to the value obtained for the 5=1 
Heisenberg model. Within a spin density wave approximation this quantity takes the value 0.93. As apparent 
from Fig. ||, the QMC data approaches smoothly this value as J/t decreases. 
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0.0 0.5 1.0 1.5 2.0 2.5 

J/t 

FIG. 6. Mean field of the two dimensional Kondo necklace model. T s (Tt) corresponds to the energy scale below which the 
bond singlets (triplets) condense. A sp denotes the spin gap. 

The quantum phase transition in the spin degrees of freedom at J c /t ~ 1.45 may be described in the framework of 
the Kondo necklace model given by: 

Hkn = t ;" r + 5/"V/') + J ' & ( 39 ) 

i,j i 

This model neglects charge fluctuations, and the spin flip processes between conduction electrons mimic the kinetic 
energy. Although the Kondo necklace model has a lower symmetry (£7(1)) than the KLM (5(7(2)) one may expect 
this model to give a reasonable description of the spin degrees of freedom at energy scales smaller than the charge 
gap. A mean-field solution is obtained in terms of bond singlet and triplet operators p2| . Both the conduction and 
impurity spins are represented by singlets, Ai, and triplets t? on the f — c bonds of the unit-cell. The bond operators 
obey bosonic commutation rules and are subject to the constraint AIAj + t^ti = 1. At the mean-field level and 
generalizing the work of Zhang et al. p3[ | to finite temperatures, one obtains the phase diagram shown in Fig. ^|. The 
condensation of singlets s — (Ai) > occurs at a temperature scale T s which, to a first approximation, tracks J. 
At J > J c the triplet excitations remain gapped and have a dispersion relation given by: uj(q) — ayl + s 2 tj(q)/a 
with a = s 2 t (l + yjl + A 2 p /t 2 s 2 ^j . Here A sp corresponds to the spin gap plotted in Fig. ||. The gap in the magnon 

spectrum at q = {tt,tt) = Q vanishes at J c /t ~ 1.4 in remarkable agreement with the QMC results. We note 
that this mean-field approach shows no phase transition in the one-dimensional case consistently with numerical 
calculations jl5|,[33]]. For J < J c the ground state has both condensation of singlets (s > 0) and of triplets at the 
antiferromagnetic wave vector ( i = V^V(i^ x ) > 0). The energy scale below which the triplet excitations condense 

is denoted by T t in Fig. ^. In terms of the KLM, the condensation of triplets (singlets) follows from the RKKY 
interaction (Kondo effect). Thus, the fact that at the mean- field level, both s and t do not vanish may be interpreted 
as coexistence of Kondo screening and antiferromagnetism in the ordered phase. We will confirm this point of view 
in the study of the charge degrees of freedom. 

B. Charge degrees of freedom 

To study the charge degrees of freedom, we compute the spectral function A(k, to) which is related to the imaginary 
time Green function via: 
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dwe"™A(fc,-cj). 



(40) 



The Maximum Entropy (ME) method is used to extract A(k,cv). Starting from the bond-singlet ground state valid 
in the strong coupling limit, one can create a quasiparticle excitation which to first order in t/J has the dispersion 
relation 



E qp (k) = 3J/4 + ty(k). 



(41) 



(7T,7r) so that the quasiparticle gap takes the value A qp = E qp (k — (jr, n)). Comparison 



2A qp > A sp in the strong coupling limit. This marks the difference to a standard band 



E q p(k) is a minimal at fc 
with Eq. @ leads to A, 

insulator which satisfies A c = A sp . In accordance with the strong coupling limit, the numerical data of Fig. |7| (a)-(c) 
show that irrespective of J/t the quasiparticle gap is defined by the k — (it, it) wave vector. Furthermore comparison 
with Fig. H shows that the inequality A c > A sp is valid for all considered coupling constants. 
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7. Single particle spectral function at T = for the ferromagnetic and antiferromagnetic KLM. We have normalized 
the peak heights to unity. The numbers on the left hand side of the figures correspond to the normalization factor. The vertical 
bars are fits to the data, see text. 



For antiferromagnetic couplings (J > 0), the spectral function shows similar features in the ordered (J < J c ) 
and disordered (J > J c ) magnetic phases thus lending support that Kondo screening which is responsible for the 
functional form of the dispersion relation at strong couplings is equally present in the ordered phase. However, 
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upon closer analysis, shadow features are seen in the antiferromagnetically ordered phase. Shadows are most easily 
understood by considering the self-energy ^2(k, uj m ) oc ^x(Q: l 'm)Go(k — q,u} m — v m ) describing electrons with 

propagator Go(k,uj m ) scattering off spin fluctuations with coupling constant g. Long-range magnetic order at wave 
vector Q — (ir, n) and staggered moment m s justify the Ansatz x(<f, v m) oc m 2 s N [38 Um fi8~ ^ for the spin susceptibility. 

The Green function is then given by: G(k,Lu m ) = 1/ ^Gq 1 (k, uj m ) — aGo(k + Q, oJ m )^j with a oc (gm s ) 2 . It is then 

easy to see that if G(k,uj) has a pole at ui then G(k + Q,uj) also has a pole at ojq, i.e. the shadow. Numerically, it 
is convenient to establish the existence of shadows by considering the imaginary time Green function. Fig. ^ plots 
(ct(r)c^) for fc = (ir, ir). At large values of rt this quantity follows an exponential law e~ rA '"'. This exponential decay 

generates the pole in A(k,uj) at u = — A qp (see Eq. (p0|). As argued above, due to the long-range antiferromagnetic 
order one expects a pole in A(k + Q, to) at lj = — A qp i.e. the shadow. As demonstrated in Fig. ||[ (ct + ^(r)cjr +( j) shows 

the same asymptotic behavior as (c^(t)cj:). Thus the low energy feature around k — (0,0) in Fig. 0c corresponds 

to the shadow of the band in the vicinity of k — (tt,tt). We note that shadow features at high energies are hard to 
resolve within the ME. Close to the phase transition in the disordered phase precursors features of the shadow bands 
are seen (see Fig. ^b). As apparent they are shifted by an energy scale which corresponds approximately to the spin 
gap. 



J /t = 1.2 




rt 



FIG. 8. ( c g +( j( r ) c fe+Q) an( i ( c g( r ) c fe) as a function of rt on a 8 x 8 lattice at T = 0. Both considered k points follow a 
e ~A qv r j aw j-gQjjj Ymes) thus confirming the existence of shadows. 

To obtain further insight into the charge degrees of freedom we will consider a mean-field theory of the KLM 
recently introduced by Zhang and Yu |54|. This mean- field theory is appealing since i) it takes into account both 
Kondo screening and magnetic ordering of the / and c electrons and ii) a phase where both Kondo screening and 
magnetic ordering emerges in a narrow region around the phase transition. 

Following Zhang and Yu |Q we write the KLM as: 

h klm ^s{k) C \ a%a + ^e(4^t - flju) (4t c ^t " i^i) + 

i 

with the constraint: f^fr^ + f\ Ji i = 1- The second term of Eq. fll^ ) describes the polarization of the conduction 
electrons by the impurity spins and leads to a magnetic instability. The third term term is nothing but a rewriting of 
the spin-flip processes: 



14 



2 ^ \ »>T i.i •'i,j.'' I .T i, 



(43) 



which are at the origin of the screening of the impurity spins by the conduction electrons. The mean- field approxi- 
mation proposed by Zhang and Yu [pij is based on the order parameters: 



{c h c n 



c£ , c? i ) = -m c e lQ ' 1 and 

1,4. 



(/t £V. . +ct f- t > = (/J O t + ct /v ,) 



-V. 



(44) 



Here Q is the antiferromagnetic wave vector, m/ and m c are respectively the staggered moments of the impurity spins 
and conduction electrons and V is the hybridization order parameter which leads to the screening of the impurity 
spins. With the above Ansatz one obtains the mean field Hamiltonian: 
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+ NJ (m f m c /4 + V 2 /2) 



(45) 



where the fc sum runs over the magnetic Brillouin zone. We note that due to particle-hole symmetry present in the 
half-filled case, the constraint of no double occupancy of the /-sites is satisfied on average: (/J./j-f + fl = 1- 



The saddle point equations, 



may then be solved [p4| 



y dmf' 



,8H 

v dm c 



(46) 



Mean Field, (n) = 2, T = 




J/t 

FIG. 9. Solution of the mean-field equations in Eq. (|46[). The solid line in (a) corresponds to the quasiparticle gap as 
obtained with the mean field order parameters plotted in (b). The dashed line corresponds to the quasiparticle gap obtained 
in the absence of magnetic ordering (i.e. we consider the solution of the mean-field equations with m c ,m/ = but V ^ 0. In 
the weak coupling limit, those solutions produce higher energy values than when magnetic ordering is allowed.) 



Solutions to the saddle point equations at T = and as a function of J/t are plotted in Fig. As apparent, 
solutions with m c) mf ^ 0, V — as well as with m c ,m/ = 0, V ^ and most interestingly with m c ,m/, V =/= 
are obtained. Each solution predicts very different functional forms for the quasiparticle dispersion relation of the 
conduction electrons. Thus by comparing with the numerical data, we can deduce which values of the mean-field 
order parameters are appropriate to best describe each phase. 
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We start by considering the spin-gap phase with J > J c . Here, magnetic order is absent and the impurity spins are 
completely screened by the conduction electrons. It is thus appropriate to set m c — rrif — but V 7^ 0. This yields 
two quasiparticle bands with dispersion relation: 

E±(k) = i (e(k) ± E{kj) , with E(k) = \j e{k) 2 + (JV) 2 . (47) 
The quasiparticle weights are given by the coherence factors: u±(k) 2 = \ ( 1 ± -^t^t )• We can use this form to fit the 



E(k) 

QMC data shown in Fig. |7] (a). As apparent, the functional form of the dispersion relation is well reproduced. 

We now consider J < J c . Here, antiferromagnetic order is present both in the conduction electrons and localized 
spins so that: m c 7^ as well as rrif 7^ 0. Following the idea that the spin degrees of freedom are frozen due to the 
magnetic ordering ordering, we set V = to obtain: 



E±(k) = ±E(k) with E(k) = \je{k) 2 + {Jm f /4) 2 (48) 

The residues of the poles of the Green function follow: u±(k) 2 = | fl ± jS^j ■ This clearly does not reproduce the 

QMC results since the very flat quasiparticle band observed numerically around k = (ir,ir) is absent (see Fig. 0(c)). 
Assuming on the other hand that magnetic ordering and Kondo screening coexist, we set V =/= to obtain four 
quasiparticle bands: 

s. 1/2 



E±,±{k) = ±^=\^E{k) ± \J E{k) 2 - ^-(m c m f /4: + V 2 ) 2 -J 2 m 2 /4:e(k) 2, j with 

E(k) = e(k) 2 + J 2 (m 2 /A + m 2 f /A + 2V 2 )/4 (49) 

An acceptable account of the numerical data is obtained by using the QMC values of the staggered moments and V 
as a fit parameter (see Fig. 0(c)). We are thus led to the interpretation that the localized spins play a dual role. On 
one hand they are partially screened by the conduction electrons. On the other hand the remnant magnetic moment 
orders due to the RKKY interaction. 

It is now interesting to consider the ferromagnetic KLM. When J < 0, Kondo screening is not present. Thus, we 
expect the appropriate mean-field solution to have m c 7^ as well as rrif 7^ but V = 0. This choice of mean-field 
parameters leads to the dispersion relation given in Eq. (p!8|). As apparent, and using rrif as a fit parameter, we can 
reproduce the QMC results (see Fig. 0d). 



IV. SPIN AND CHARGE DEGREES OF FREEDOM AT FINITE TEMPERATURE 



The aim of this section is to define relevant energy scales for both spin and charge degrees of freedom as a function of 
J/t. In doing so, we will discuss the behavior of the optical conductivity, staggered spin susceptibility, single particle 
spectral functions as well as specific heat as a function of temperature. We will put the emphasis on the behavior of 
those quantities at the spin and charge energy scales. 



A. Spin and Charge energy scales. 

To define the charge scale, we consider the charge susceptibility % c = -j% ((N 2 ) — (-/V) 2 ) where N corresponds to 
the particle number operator. It suffices to consider only the conduction electrons since the /-electrons are localized 
and have no charge fluctuations. Since we are discussing the temperature dependence of Xc, let us recall the high- 
temperature result: 

^('-sM^ 2 ))' (50) 

From that behavior, it appears that J\ c will exhibit some approximative scaling form as a function of T / J only for 
large J/t » 8/V3 ~ 4.62. In Fig. [To[ we find consistent results at high temperature with (|50|). 
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FIG. 10. Charge susceptibility J\c vs T/J for various couplings on the L — 6 lattice. For very large J, we obtain good 
agreement with the large- J expression (^l|) plotted with a dashed line. 

We can define a characteristic charge temperature, Tc, in a precise way by looking at the maximum of Xc- In the 
weak coupling limit, our numerical results are consistent with Tc ~ J (See Figs. |l0| and [l9]). In the large J limit (J > 
bandwidth), the physics becomes local and one can consider decoupled sites. For each site, there are only 8 states to 
take into account for computing the grand-canonical partition function: the singlet state, the three-fold degenerate 
triplet, the four-fold degenerate 5=1/2 state containing cither an empty conduction site or a doubly-occupied one 
and with the two different spin configurations. In this limit, A sp = J, A qp / J — 3/4 and 

4 

Xc= ^4 + 3e-^/ 4 + e 3 ^/ 4 (51) 

which exhibits a peak at Tc = 0.386J. Hence, and apart from different numerical prefactors at weak and strong 
couplings, Tc scales as J in both limits (see Fig. |Ts| ). 

To best understand the meaning of the charge scale, we consider the real part of the optical conductivity as obtained 
from the Kubo formula, a(u>,T). This quantity is related to the imaginary time current-current correlation functions 
via: 

(J(r)J(O)) = / duK(w,r)a(u;,T) with K(w,t) = ~ 1 _ e -p u - (52) 

Here J is the current operator along the x or y lattice direction and () represents an average over the finite-temperature 
ensemble. The above inverse Laplace transform, to obtain the optical conductivity is carried out with the ME |is|| 
method. The default model is chosen as follows. We start at high temperature with a flat default and then, for lower 
temperatures, we take as default the result obtained at the temperature just above pq| . This allows us to obtain 
smoother results but emphasizes the fact that the ME method depends on the default which is used. 

The overall features of the conductivity are shown in Fig. ^ for a given J. At high temperatures, there is only a 
very broad lorentzian Drudc peak. By lowering the temperature, we first observe an enhancement of the Drude weight 
as expected for a metal. At temperatures scales lower than Tc, there is a transfer of spectral weight from the Drude 
peak to finite frequencies and finally, at very low temperatures, we observe the opening of an optical gap related to 
the quasi-particle gap observed in the single particle density of states (DOS). 
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FIG. 11. Optical conductivity a(ui,T) versus w on a logarithmic scale for J/t = 0.8 and various temperatures (At = 0.2, 
L = 8 lattice). The peak height has been normalized to unity and the normalization factor is listed on the left hand side of 
the figure. As the temperature is decreased below the charge scale, Tc/t ~ 0.16, spectral weight is transferred from the Drude 
peak to finite frequencies. 



The resistivity is defined as p(T) = 1/ct(0,T). In Fig. 
located at approximately Tc- Thus, we will conclude that 
electrons is enhanced while decreasing temperature due to 



|l2| , we plot p(T) for various J. We observe a minimum 
Tq corresponds to an energy scale where scattering of the 
the screening of magnetic impurities. 




FIG. 12. Normalized dc resistivity p as a function of T/J for various couplings. We have checked for J/t — 0.8 that the 
results do not depend on At. To a first approximation, and taking into account the scatter of the data at J/t = 1.6, the 
temperature of the minimum in p tracks Tc- Tc/J ~ 0.2, 0.2, 0.25 for J/t = 0.4, 0.8 and 1.6 respectively. 
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FIG. 13. Normalized local spin-spin correlation function C = (5/ ■ S c )(T)/\{Sf ■ S C )(T — 0)| as a function of T/J for various 
couplings. To a first approximation, the temperature scale of the onset of correlations tracks J . 

This scenario is reinforced by the behavior of the local spin-spin correlation function C = (Sf ■ S c ) (T) /\(Sf S c ) (T = 
0)| plotted in Fig. [l3|. As the temperature is lowered, this quantity decreases indicating the formation of local singlets. 
Since the curves are almost identical as a function of Tf J and for various couplings, we deduce that the typical energy 
scale is J and that the formation of those singlets are responsible for the enhancement of the resistivity which occurs 
at a similar temperature. 

Before considering the characteristic energy scale for the spin degrees of freedom we comment on the relation 
between the optical gap - as obtained from the low temperature conductivity data - and quasiparticle gap (see Sec. 
III). They are not directly related since optical transitions involve only zero momentum transfer. Starting from the 
hybridization picture, we can represent the band structure as in Fig. |l4]. Generalizing this figure to 2D, we clearly 
see that the smallest optical gap is at k — (it, 0) (or equivalent points) and is larger than the charge gap A c ~ 2A qp . 
We recall that the quasi-particle gap corresponds to a transfer from a particle at k — (tt,tt) in the lower band to the 
chemical potential. We then expect from that naive argument A op t > Ac — 2A gp . More precisely, we can relate the 
optical gap to the gap at k = (0, n) as observed in Fig. U3. 



E(k) 
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FIG. 14. Schematic ID band structure showing the hybridized bands, k varies for k = to k = -k and as apparent charge 
gap Ac ~ 2A 9P is smaller than the optical gap A op t. 
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FIG. 15. Various gaps as a function of J. We have considered low enough temperatures so as to reproduce ground state 
results and L = 6. The quasi-particle gap is obtained from the DOS, the gaps at (jt,tt) and (0,0) are seen in the spectral 
functions at those points, the optical gap, A CT stems from the optical conductivity and finally, the spin gap A sp is taken from 
Ref. As apparent, in the weak coupling limit, where the quasiparticle dispersion is very flat along the k = (7r,7r) to 

k — (0, 7r) direction, the optical and charge gaps are comparable. (See Fig. [j]) 



To define a characteristic energy for the spin degrees of freedom, we compute the uniform spin susceptibility, 
Xs — j2 (( m z) — ( m z) 2 )- Here, m z = ^2i(n^^ — 7ijjJ with a — ct ^ a + ft /| a . In order to observe magnetic 
properties, it can be necessary to go to very low temperatures when J is small. With our algorithm which is free from 
the sign problem, we can go down to T = O.Oli for L = 6 or T = 0.02£ for L = 8. 
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FIG. 16. Uniform spin susceptibility J\ s as a function of T/J for various couplings and lattice sizes. The high-temperature 
expansion (|53[) is shown with filled circles and the large-coupling expansion (154) is plotted in dashed line. 



In a high-temperature expansion, \s takes the form: 



= — (\ - — 

Xs 8T V 6T 



(53) 



From this expansion, one expects to observe a scaling property Jx s — f(T/J) in this regime. This is indeed what 
is observed in Fig. |l6| for T/J> 0.6. We define the magnetic characteristic temperature Tg via the position of the 
maximum in Xs- At large coupling, the physic of the Kondo lattice becomes local. In that limit, the susceptibility is 
easily computed (see a similar calculation for the charge susceptibility in Eq. (51)) and takes the form : 



Xs 



1 + 2e-' 3 - / /4 



4 + 3 e -PJ/i + e 3pj/i ' 



(54) 



which exhibits a maximum at Ts — 0.453 J. In contrast, for smaller J, the position of the maximum clearly increases 
more slowly than J (see Fig. |l7|). As apparent from Fig. 17 and for the considered values of J/i, Ts scales approximately 
as J 2 . 



J 2 Xs 6 - 
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FIG. 17. Uniform spin susceptibility J 2 Xs as a function of T/J 2 for various couplings and L — 8. For J < 6, the maximum 
Ts/ J 2 is roughly constant. 

Fig. |l9| plots Ts as a function of J. In the large coupling region, we have excellent agreement with the expansion 
of Eq. Q54|); but, for couplings up to ~ 5 (or bandwidth which is the physical case), Ts is well fitted by ~ J 2 . 

The meaning of the energy scale Ts is elucidated by considering the spin susceptibility at the antiferromagnetic wave 
vector Q = (jr, it). This quantity measures the antiferromagnetic correlation length. Indeed, writing the spin-spin 
correlation functions in space and imaginary time as S(f,r) = Ae~xjp(ir ■ Q) exp(— r/£) exp(— r/£ T ), we find that the 
staggered susceptibility Xs(Q) = Xf dr J drexp(—iQ ■ f)S(r, r) ~ £ D £, T in D dimensions. For the Heisenberg model |, 
the dynamical exponent z defined by £ r ~ £ z is equal to 1 plJ3^]. We then obtain in our case Xs(Q) ~ £ 3 - 



2 Since the charge degrees of freedom are gapped, we expect that our model is in the same universality class as the 0(3) model. 
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FIG. 18. Staggered spin susceptibility Xs(Q) f° r various couplings and sizes. Since Xs(Q) ~ £ 3 > we can extract the behavior 
of the antiferromagnetic correlation length f . T s ~ 0.017, 0.05, 0.22 for J/t = 0.4, 0.8, 1.6. 

Xs(Q) is plotted in Fig. As apparent and for the considered J/t range, the energy scale Ts marks the onset 
of s/iort-range antiferromagnetic fluctuations. At low temperatures in the ordered phase, one expects £ to grow 
exponentially as a function of decreasing temperature. On the other hand, in the spin gap phase, J/t > 1.45, the 
antiferromagnetic correlation length saturates to a constant |37|] . 
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FIG. 19. Characteristic spin Ts and charge To temperatures as defined by the maximum of 
for L — 6 and L = 8 when available. At large J, the asymptotic behavior of T$ is 0.453J (Eq. (p 
adjustable parameter; at small J, Ts is well fitted by the form J 2 /II (dashed line). 
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The results of this section are summarized in Fig. We have defined both a charge, T c , and a spin, Tg, scales. 
The charge scale corresponds to the onset of enhanced scattering as a function of decreasing temperature due to the 
screening of the magnetic impurities. Apart from different numerical prefactors Tc scales as J in the weak and strong 
coupling limits. From weak to intermediate couplings the spin scale defines the onset of short range antiferromagnetic 
fluctuations and follows a J 2 law in agreement with the energy scale associated to the RKKY interaction. At strong 
couplings, Tg tracks the spin gap. We note that we find good agreement with exact diagonalizations studies at finite 
temperatures p8j . This approach is however limited to very small cluster sizes and consequently to high temperatures 
and/or large values of J/t where the local approximation becomes valid. Thus those studies cannot extract the 
behavior of Tg in the weak coupling limit. 



B. Temperature dependence of spectral functions and origin of quasiparticle gap. 

The origin of the quasiparticle gap in the strong coupling limit is the formation of Kondo singlets. In the weak 
coupling limit the situation is not a priori clear. In the mean-field theory presented in Fig.|9| and retaining only 
Kondo screening, we obtain an exponentially small gap corresponding to the dashed line in Fig K3a. On the other hand 
retaining only magnetic ordering, the quasiparticle gap takes the value J/At in good agreement with the numerical 
data. We note that an exponentially small gap is equally obtained with (i) Gutzwiller approximation p2| ; (ii) 
dynamical mean-field theory [jioj; (iii) 1/N expansion ]39[ | since those approximations neglect magnetic fluctuations. 
In this section, we argue that at or slightly below Tc a small gap emerges leading to the quasiparticle dispersion 

relation | ^e(^) ± \j e 2 (k) + A 2 ^ /2 and that the quasiparticle gap of order J is formed only at Tg. 

We start by considering the integrated DOS, N(u>) obtained with the ME method. Results are shown in Fig. at 
J/t = 0.8. In the vicinity of the charge scale, Tc — O.lQt, one observes a reduction of spectral weight at the Fermi 
energy. Within the mean-field approximation of the KLM presented in Eq. (^) and (p4|), this dip in the DOS of the 
conduction electrons follows directly from the occurrence of Kondo screening, i.e. V ^ 0. Hence this feature in N(u>) 
at Tc stands in agreement with our interpretation of the charge scale Tc- As the temperature is lowered below Tc, 
the density of states at the Fermi level is further depleted and a gap opens in the low temperature limit. 
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FIG. 20. Integrated DOS for J/t — 0.8 on the L = 8 lattice for various temperatures shown on the plot. The peak height is 
normalized to unity and normalization factor is listed on the left hand side of the figure. 

In order to gain more insight into the distribution of spectral weight, it is convenient to compute the momentum- 
dependent DOS, A(k,u>). The integrated density of states merely corresponds to the sum over all k of A(k,u>). 
Therefore, we expect the same behavior by decreasing the temperature but we have more information on the dispersion 
relations of the excitations for example. Fig. |2l| a plots A(k, lo) again for J/t = 0.8 and at a temperature T = 0.083t 
corresponding to T$ < T < Tq- For comparison, we have included the T = data (see Fig. pT)b) As apparent the 
substantial spectral weight of the k points on the non-interacting Fermi line i.e. k — (0, 7r), (it/2, tt/2) has shifted to 
lower energies. This is the origin of the decrease in spectral weight observed at the Fermi level in the integrated DOS 
at T ~ Tc- However, the flat dispersion relation around k = (tt, 7t) - with significantly less spectral weight - remains 
pinned at the Fermi level. The dominant features of the quasiparticle dispersion relation are well reproduced by the 

fit: \e{%) ± y / e 2 (fc) + A 2 ^j /2 with A = QM. This value of A produces a quasiparticle gap A qp = A 2 /16t ~ 0.016t 

which lies beyond our resolution. As seen in Fig. [2l|b, A qp — 0.28 ± 0.02 in the zero temperature limit. 
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FIG. 21. (a) Spectral functions for J/t = 0.8 and (3t — 12 (At = 0.2, L = 8 lattice). Normalization factors are written on 
the vertical axis. At this temperature, Ts < T < Tc = Q.2t, there is a clear formation of hybridized bands with quasiparticle 
gap lying beyond our resolution. The vertical bars correspond to a fit of the data (see text). For comparison we have included 
the T = results (b). 

Since the quasiparticle gap is determined by the k = (jr, tt) wave vector we concentrate on this fc-point to analyse the 
temperature evolution. As apparent in Fig. |22| at J/t = 0.8 the quasiparticle gap of order J is formed approximately 
at the magnetic scale Ts — 0.05t. Since the model is particle-hole symmetric A(k,oj) = A(k + Q, —w). Thus the fact 
that the peak splits symmetrically around the Fermi energy confirms the presence of shadow bands. In the spin gap 
phase the quasiparticle gap originates solely from Kondo screening. In the mean-field approximation presented in 
Eq. ( J4^ ) and (Q) and retaining only Kondo screening, the quasiparticle gap will grow continuously as a function of 
decreasing temperatures below the charge scale. This merely reflects the temperature dependence of the mean-field 
order parameter V. Precisely this behavior is seen in Fig. ^2] at J/t = 1.6. 

J/t=0.8 J/t=1.6 




J/t = 0.8, the quasiparticle gap of order J/4t opens at a temperature scale comparable to T$ = 0.05. In the spin gap phase, 
at J/t = 1.6 the quasiparticle gap grows smoothly as a function of decreasing temperature. 

The evolution of the quasiparticle gap as a function of temperature is equally seen in the charge susceptibility. At 
low temperatures one expects \c = /3exp(— A qp /3). As apparent from Fig. [23l it is only below Ts that the data follows 
the above exponential form. 



J/t = 0.8, L = 8 
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FIG. 23. Xc/fl a t J/t = 0.8. The solid line corresponds to exp(— A 9P /3) where quasiparticle corresponds to the value 
obtained with T = simulations (See Fig. ^). Only below the spin scale, Ts, do the data follow the above exponential form. 
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C. Specific heat 

Finally, we consider the specific heat which contains information on both spin and charge degrees of freedom. In 
principle one can obtain the specific heat C V (T) by direct calculation of the fluctuations of the internal energy E(T): 
C V (T) = jjdE/dT = -^((H 2 ) - (H) 2 )/T 2 . However, this method produces very poor results at low temperatures. We 
have thus used a ME method to compute C v as proposed in ]io| . In Fig. ^J, we show C V (T) as well as the uniform 
spin an charge susceptibilities for various couplings as a function of temperature. 

At J/t = 0, the specific heat is given by the sum of a delta function at T — for the localized spins and the specific 
heat of free fermions. By switching on the coupling, they are combined to form a two-peak structure. The broad peak 
at high temperature T ~ t is almost independent of the coupling J and is rather similar to the free electron gas. The 
sharp peak at lower temperatures strongly depends on the exchange constant. It shifts toward higher temperatures 
and becomes broader with increasing J/t. The location of this peak tracks the magnetic scale T$ indicating that its 
origin comes from the spin excitations. In the spin gapped phase, we note that the overall features of C v agree with 
the ID case Q. 
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V. SUMMARY AND CONCLUSIONS 



We have presented a detailed numerical study of ground state and thermodynamic properties of the ferromagnetic 
and antiferromagnetic half-filled KLM model on a square lattice. From the technical point of view, we have described 
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and used an efficient (i.e. free of the minus-sign problem) auxiliary field QMC method to investigate the model. Both 
finite and ground-state algorithms were discussed. The approach is by no means restricted to the KLM and may be 
applied to investigate models such as the half-filled two channel Kondo lattice or various forms of depleted Kondo 
lattices in which the impurity spins are removed in a regular or random way. However, we are tied to particle- hole 
symmetry since only in this case can we avoid the minus-sign problem. 

In two dimensions the KLM shows a quantum phase transition between antiferromagnetically ordered and disordered 
states. This transition occurs at J/t = 1.45 ± 0.05. The magnon dispersion evolves smoothly from its strong 
coupling form with spin gap at Q = (tt,tt) to a spin-wave form in the ordered phase. The transition may be well 
understood in the framework of a bond-operator mean- field approximation of the Kondo necklace model ||33|] . Here, 
the disordered phase is represented by a condensation of singlets with an energy gap at Q = (7r, 7r) for magnon 
excitations. At and below the critical point the spin gap closes leading to a condensation of both singlets and triplets 
at the antiferromagnetic wave vector. The system remains insulating. To a first approximation and as in the one- 
dimensional case, the quasiparticle gap scales a | J| irrespective of the sign of J. In contrast, the quasiparticle dispersion 
relation shows marked differences between ferromagnetic and antiferromagnetic couplings. For antiferromagnetic 
couplings the quasiparticle dispersion always has a structure which follows the functional form obtained in the non- 



interacting PAM: E±(k) = i ^e(^) ± y s(k) 2 + A 2 J. This functional form is obtained in various approximations 

|54],[l2| which take into account Kondo screening but neglect magnetic ordering. In the antiferromagnetic phase the 
above dispersion relation is merely supplemented by shadow features. One obtains a four band structure which is well 
reproduced by mean-field theories which produce phases with coexistence of magnetic ordering and Kondo screening 
P4f . Thus, in the ordered phase screening of the impurities is not complete. The remnant magnetic moments order 
due to the RKKY interaction. Although we cannot dope the system - due to severe minus-sign problems - it is 
tempting to assume a rigid band picture and to describe the doped phase by shifting the chemical potential into 
the conduction band. Since the quasiparticle gap is determined by the k = (±7r,±7r) points, the Fermi line will 
consist of hole pockets around those points and one expects the Luttinger volume to account both for localized and 
conduction electrons. Furthermore, since the band is very flat around those points a larger effective mass is anticipated. 
Ferromagnetic couplings show a different behavior. In this case, Kondo screening is absent but the RKKY interaction 



present. The quasiparticle dispersion is well fitted by the form E±(k) = ±y e(k) 2 + A 2 corresponding to free electrons 
subject to an external staggered magnetic field. In this case, again assuming a rigid band picture, doping produces 
a Luttinger volume containing only the conduction electrons. This contrasting behavior of the Luttinger volume for 
the ferromagnetic and antiferromagnetic KLM is reproduced in the limit of large dimensions [[llj . 

From the finite temperature simulations, we can define spin, T5, and charge, Tc, energy scales by locating the 
maximum in the charge and spin susceptibilities. From weak to intermediate couplings the spin scale follows a J 2 law 
in agreement with the energy scale associated with the RKKY interaction. At strong couplings Tg cx J. In contrast 
both in the weak and strong coupling limit Tc ~ J- In the range where Ts oc J 2 , the staggered susceptibility shows a 
marked increase at T ~ T<j. Hence, in this range T$ corresponds to the onset of antiferromagnetic fluctuations. On the 
other hand, the charge scale determines to a first approximation the minimum in the resistivity. Furthermore, at Tc 
antiferromagnetic intracell correlations between the /- and c- electrons are formed and a dip in the integrated density 
of states, N(oj), at the Fermi level is observed. Thus, this scale marks the onset of enhanced scattering originating 
from the screening of the magnetic impurities. In the limit of infinite dimensions, a similar behavior in the charge 



degrees of freedom is seen but at a much smaller energy scale, T - e -i/zJN(e s ) ^ Apart form a factor 1/2 in the 
exponent To corresponds to Kondo temperature of the single impurity problem. In one dimension, a dip in N(tu) is 
observed at an energy scale larger than the spin gap which scales as e~ 1 / aJN ^ es ^ in the weak coupling limit (with a 
numerical estimation ofl<a<5/4 jl5| or a = 1.4 [fL6|). 

In the weak coupling limit, one can analyze the single particle spectral function at various temperatures. Our results 
show that the quasiparticle gap of order J is formed only at the magnetic energy scale. Thus one can only conclude 
that the quasiparticle gap at weak couplings is of magnetic origin. In contrast at strong coupling, the quasiparticle 
gap originates from Kondo screening. The above stands in agreement with arguments and numerical results presented 
for the one dimensional case ]15| , |l6| ]. At weak couplings in ID, the spin gap becomes exponentially small. Hence, 
the time scale associated with magnetic fluctuations is exponentially larger than the time scale relevant for charge 
fluctuations which is set by t. The conduction electrons thus effectively feel a static magnetic ordering. In ID and in 
2D in the presence of particle- hole symmetry, nesting of the non-interacting Fermi surface is present. At a mean- field 
level and in the presence of magnetic ordering, this leads to a quasiparticle gap A qp — J/4. In 2D one may alter the 
shape of the non-interacting Fermi surface so as to avoid nesting by introducing a small nearest neighbor hopping 
matrix element. In this case, the mean-field approximation does not produce an insulating state in the presence 
of antiferromagnetic ordering. Since nesting is related to particle-hole symmetry, we cannot address this question 
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in the QMC approach due to severe sign problems. Hence it is worth paying particular attention to our results at 
weak couplings and Tc > T > T$ before antiferromagnetic correlations set in. Here, Kondo screening is present but 
antiferromagnetic correlations absent. In this temperature range, A(k,uj) shows a dispersion relation following that 
of hybridized bands with quasiparticle gap lying beyond our resolution. 

We have equally computed the specific heat, C v . This quantity shows a two-peak structure. The broad high energy 
(T ~ f) feature stems from the conduction electrons. The low energy peak is very sharp in the ordered phase and 
tracks Ts- It is hence of magnetic origin. 

Finally we discuss the relationship of our results to experiments. Let us first concentrate on Ce3Bi4Pt3. At 
T = 100 K the effective magnetic moment of Ce ions starts decreasing E3]. At higher temperatures the Ce ion 
has a next to fully developed moment (i.e. J = 5/2 as appropriate for Ce 3 ^). At the same temperature scale the 
real part of the optical conductivity shows a reduction of spectral weight in a frequency range of 39meV or 450 K 
Those results imply that the opening of a gap is related to the screening of magnetic impurities and hence, 
the KLM seems to be an adequate prototype model for the description of this class of materials. The above described 
temperature evolution is precisely seen in our numerical simulations. At T ~ Tc and at weak couplings, the optical 
conductivity shows a transfer of spectral weight from low frequencies to frequencies well above Tc (Fig. |ll"| ) . Screening 
of the magnetic moments start equally at T ~ Tc (Fig. |l3|). For the above material, the optical gap is estimated by 
A CT = 39meV Q and photoemission experiments suggest a quasiparticle gap A qp = 20meV Q. At a temperature 
scale T ~ 25 K a gap in the magnetic excitation of A sp = 12meV is observed. Those small energy scales imply 
that small values of J/t should be considered. The gaps equally satisfy the relation A c ~ 2A ?p > A sp as obtained 
in the KLM. Hence, one should place this material in the parameter range J > J c which in our calculations seem 
rather large in comparison to the small charge gap observed in experiments. However, one should keep in mind that 
J c may be sensitive to the properties of the non-interacting Fermi surface. In particular nesting - which is present 
in our calculation - will certainly enhance the value of J c . We now turn our attention to CeNiSn. CeNiSn has a 
transport gap roughly an order of magnitude smaller than Ce3Bi4Pt3, and hence - assuming a KLM description of the 
material - should correspond to smaller values of J/t in comparison to Ce3Bi4Pt3. This smaller value of J/t leads 
to signs of magnetism. Indeed along the a-axis of the orthorhombic structure, CeNiSn shows a peak in the magnetic 
susceptibility at 12K. At the same energy scale, an anomaly is seen in the specific heat p8| . This seems consistent 
with our results. 
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